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Practical  and  theoretical  considerations  in  computing  first  passage  time 
statistics  are  considered.  We  are  motivated  by  first  passage  times  as  models  of 
failure  times. 


In  particular,  let  X(t)  be  a  diffusion  on  [0, r ]  with  a  reflecting  boundary  at  0. 
Denote  by  tt  the  time  of  first  passage  to  level  r,  ie,  rr  =  inf  {t  >  0  ,  X(t)  >  rj, 
and  let  w(x,t )  be  its  tail  probability  function  conditional  on  X(0)  *=  0  ,  ie 
w(x,t)  =»  P{  Tr  >  t  |  JST(O)  =  x}  =  Px{rr  >  f}  . 


In  Section  1,  the  relevance  of  first  passage  time  distributions  as  failure  time 
models  is  indicated/(ef  [8}).'  Also,  the  spectral  series  expansion  solution  to  the 
backward  equation  is  introduced. 


In  Section  2,  algorithms  for  approximating  w(x,t)  are  obtained.  In 
particular,  the  infinite  spectral  expansion  for  w(x,t)  is  approximated  by  an  n- 
term  sub-expansion  which  matches  the  first  n—  1  moments.  Proofs  validating 
the  spectral  expansion  and  the  related  approximation  scheme  are  given  in  the 
Appendix. 


In  Sections  3  and  4,  methods  are  given  for  obtaining  the  eigenvalues  and 
first  passage  moments,  necessary  for  computing  approximations  to  w(x,t)  .  In 
Section  5.  computational  issues  related  to  calculating  the  moment  generating 
function  are  considered. 


Sections  6  and  7  include  theoretical  complements  about  first  passage  times. 
In  particular,  the  moment  generating  function  is  shown  to  possess  an  interesting 
representation  having  exponential  form  (cf  equation  (7.1)  ).  This  exponential 
representation  is  related  to  asymptotic  expansions  used  in  analyzing 
perturbations  of  certain  second-order  differential  equations. 
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1.  Introduction  And  Motivation 


Let  {fi  ,  F  ,  P}  be  a  probability  space.  Suppose  that  Ft  is  a  filtration  of  F 
and  that  (VF(i),  t  >  0}  is  standard  Brownian  motion  adapted  to  Ft. 

Let 

A  =  ^p~D2  +  n(x)D  (1.1) 

be  the  infinitesimal  generator  of  a  diffusion  X(t)  on  [0  ,  r]  satisfying 
dX(t)  =  a{x)dW{t)  +  n(x)dt  , 

with  a  reflecting  boundary  at  0  and  absorption  at  r  <  oo  ,  and  where  cr^x)  >  0 
and  n{x)  are  continuously  differentiable  on  [0  ,  r]  . 

Define  the  stopping  time  rr  by 

rr  =  inf  {s:  $>0  ,  X(s)  >  r} 

and  the  moment  generating  function  %(x,y)  by 

^(x,y)  =  E{etr>  |X(0)  =  x}  =  Ez[e'r’} 


1.1  First  Passage  Times  As  Failure  Times 

Our  motivation  for  studying  first  passage  time  distributions  is  their 
relevance  to  modeling  of  failure  times.  Indeed,  this  paper  continues  the  line  of 
development  initiated  in  [8],  where  a  stochastic  process  is  used  to  model  system 
state,  ie,  wear-and-tear,  and  failure  occurs  when  either  a  traumatic  killing  event 
occurs  (killing  events  happen  with  rate  £(x)  in  state  x),  or  the  system  is  retired 
when  wear-and-tear  reaches  some  predefined  threshold  (ie,  a  first  passage 
occurs). 

For  example,  if  system  state  is  modeled  as  Brownian  motion  with  positive 
drift,  then  first,  passage  to  a  specified  threshhold  has  an  inverse  Gaussian 
distribution.  This  first  passage  distribution  has  been  successfully  applied  to 
numerous  problems  to  obtain  good  fits,  (cf  Jorgensen  [4]  ). 

A  related  but  parallel  line  of  development  is  explored  in  Wenocur  [ll]  , 
where  the  killing  time  distribution  of  Brownian  motion  with  quadratic  killing 
rate  is  calculated.  ° 

Our  aim  in  this  paper  is  to  study  first  passage  time  distributions,  where  the 
system  state  process  is  a  general  diffusion  with  reflection  at  the  origin  and 
absorption  at  r  <  oo.  That  is,  the  system  state  evolves  as  a  diffusion,  and 


□ 

□ 


ty  Codes 
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failure  occurs  at  the  epoch  of  first  passage  (or  absorption)  to  level  r  . 

In  future  work,  we  intend  to  explore  the  practical  ramifications  of  employing 
the  computational  methods  suggested  here  to  evaluate  interesting  first  passage 
times  statistics. 


1.2  Backward  Equation  For  First  Passage  Time  Distribution 

Let  w(x,t)  denote  the  tail  of  the  first  passage  time  distribution,  ie, 
w(x,t )  =  Px{tt  >  t}  . 

The  backward  differential  equation  for  w(x,t)  is 

ishJl  »  Jstedl  +  =  Av(Xit)  (1.2) 

for  (x,t)  e  (0 ,r)  X  (0,oo)  ,  with  boundary  conditions 

w(x,0)  **  1  for  0  <  x  <  r  ,  and  for  all  t  >  0  w(r,t )  =0  and  —  q  . 

ox 

For  a  derivation  of  this  equation  and  other  related  quantities  see  [5],  pp  222- 
224. 

1.3  The  Spectral  Representation  for  w(x,t) 

The  follov'ing  representation  for  w(x,t)  is  valid  whenever  cr( x)  and  p(x)  are 
sufficiently  smooth. 

w(x,t)  =  £  ck  e~akt<j>k(x)  (1.3) 

x-i 

where  ak  and  <f)k  are  eigenvalues  and  eigenfunctions,  and  ck  are  generalized 
Fourier  coefficients,  all  defined  below  (This  representation  is  proved  in  Section 
8.1). 

The  (  <bk  ,  k  >  1}  are  eigenfunctions  of  A  corresponding  to  the  eigenvalues 
H  ,  *  >  l},  ie, 

A<t>k  =  — ^ ock<j>k  , 


ck  “  /  4>k{x)p{x)dx  , 


JUVS? 

'  S.t-  •  * 


■ 


where  p(x)  is  given  by 
p(x)  -  2ir{x)/o2{x) 
and  n(x)  is  given  by 


(1.4) 


tt(x)  =*  exp/  2p(u)fo1(u)du  .  (1-5) 

o 

In  general  an  arbitrary  function  /  G  L2{p)  will  have  a  Fourier  type 
expansion,  ie, 

/  *  S  ck<t>k 

*-i 


where  equality  is  interpreted  in  the  L2(p)  sense  and 

r 

C*  *  / / {x)M*)p{x)dx 
o 

Remark:  In  the  sequel  it  is  assumed  that  A ’s  eigenvalues  form  a  complete  set  in 
L2(p).  The  completeness  of  A’s  eigenfunctions  can  be  assured  by  certain 
regularity  conditions  on  the  infinitesimal  parameters  cr(x)  and  p(x)  .  For 
example,  cr(x )  >  0  and  the  continuity  of  <r"(x )  and  p’{x)  are  sufficient 
conditions.  See  [10,  chap  l]  for  more  details. 


1.4  A  Generalization 

This  paper  is  primarily  concerned  with  computing  first  passage  time 
statistics.  In  [8],  as  alluded  to  in  (1.1),  a  general  reliability  model  was  proposed 
in  which  system  failures  occur  when  either  system  wear-and-tear  reaches  some 
maximum  permissible  level  (ie,  a  first  passage  occurs)  ,  or  when  some  killing 
event  happens  (such  killing  events  occur  with  rate  k{x)  in  state  x).  Under  this 
model  w(x.t)  satisfies  the  following  equation: 


,  k(x)w(xmx)S^Ai  +  ih)  =  BwM , 

with  the  same  boundary  conditions  as  (1.2)  ,  and  where 

Bf{x)  -  ~<T 2(x)/"(x)  +  p{x)J'{x)  +  k{x)J{x)  . 

It  is  possible  to  solve  for  w(x,t)  and  related  quantities  with  methods  very 
similar  to  those  presented  here. 
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2.  Approximating  The  First  Passage  Time  Distribution 


The  above  discussion  might  suggest  that  solving  for  w(x,t )  is  pretty 
straightforward.  But  generally  eigenvalues  and  especially  eigenfunctions  are 
difficult  to  obtain.  However  the  problem  of  approximating  w(x,t)  can  be 
approached  by  the  method  of  moments.  One  technique  is  to  calculate  the  first 
three  moments,  and  then  use  the  Pearson  curve  fitting  method  (cf  [9]).  This 
method  is  computationally  feasible,  and  the  Pearson  family  of  curves  includes 
some  important  first  passage  distributions,  such  as  the  gamma  distribution  (cf 
[l]).  The  merits  of  this  approach  will  be  studied  in  a  forthcoming  paper. 


2.1  Given  n  Eigenvalues  And  n— 1  First  Passage  Moments 

A  more  computationally  intensive  approach,  but  one  founded  on  stronger 
theoretical  grounds,  is  the  following.  Suppose  that  n  moments 
{Mk(x,r)  ,  1  <  k  <  n  }  are  known,  where  Mk(x,r)  =>  E1  [  r*  ],  as  well  as  the 
first  n  eigenvalues  {a*  ,  1  <  &  <  n  }  .  Then  use  a  finite  sum  in  place  of  the 
infinite  sum  in  equation  (1.3).  In  particular,  approximate  w(x,t)  by 

w„(*.0  -  £  pf\*) 

j- 1 

where  p1"’  s  (p[n)  ,  .  .  .  ,p^})  satisfies  for  0  <  k  <  n—  1 
*  ,  Mk[x.r) 

ZJ  P]  (x)  (Vj)k  - - -  Where  Hj  =  1  /<*j  (2.1) 

y-i  K- 

Ideally  we  want  wn(x,t)  to  be  a  distribution  function,  ie,  wn{x,t)>0  and 
wn(x,s+0  >  wn(x,t)  whenever  s  >  0.  It  is  not  clear  that  solving  (2.1)  always 
produces  such  a  function.  This  issue  requires  further  investigation. 

2.1.1  Obtaining  The  Weighting  Factors 


The  weighting  factors  {p*("’  ,  1  <  k  <  n  }  in  (2.1)  above  are  obtained  by 
solving  the  following  linear  system: 


1 

Ml 

0 

K 

1 

m2 

m2 

1 

...  M.V 

■  ■  ■  M.v 

Pr 

pr 

H 

i 

m  j 

m2 

X,"-' 

Mr1 

.  .  .  ,,n-l 

M/v 

mn-l 

where  Hj  *  l/or;-  and  mk  *  Mk[x,r)/k\ 


Observe  that  the  matrix  {fij  ,  1  <  j  <  n  ,  0  <  ,?  <  n}  is  none  other  than  the 
transpose  of  the  celebrated  Vandermonde  matrix.  Cramer’s  equation  gives  the 
following  formula  for  p£%)  . 

Pk’]  =  E  n  (M>  -  Mi)  (2-3) 

i  <  ;  <  #  !<;'<» 

;*F 

where  gr  are  the  signed  symmetric  functions  defined  as  follows: 

0o(ai>a2>  ‘  '  »am)  =  1 

and  for  r  >1 

gr(ava2,  ■  ■  ■  ,am)  -  2  a,- a,-2  •  •  • 

1  <  *  i  <  *  a  <  •  <  i,  <  m 

It  is  shown  in  Section  (8.2)  that  p^*’  —  pk  —  0(n~2)  . 

2.2  Given  2n—l  Moments  Only 

Suppose  that  the  first  2n—l  moments  have  been  determined.  It  is  possible  to 
approximately  determine  the  first  n  eigenvalues  by  solving  the  following  system 
of  equations  for  At,  ,  1  <  i  <  n. 


1  1  1 

r 

1 

Ml  M  2  ...  Mn 

2  o  o 

Ml  Mo  •  •  •  Mn 

pV 

. 

mj 

m2 

2n-l  ..  2n  — 1  .  .  .  ,.2n-l 

Ml  Mo  Mn 

where  fXj  >  0  . 

Mn’’ 

m2n  —  l 

One  approach  is  to  solve  for  p(’’1  in  terms  of  {pk  ,  l  <  k  <  n}  and 
(m0  .  ■  .mn_!}  as  detailed  above,  and  then  use  the  remaining  n 

constraints  to  determine  the  {pk  ,  1  <  k  <  n}.  This  reduced  system  can  then  be 
solved  using  mathematical  programming  techniques,  eg,  approximate  Newton- 
Raphson  techniques.  The  numerical  stability  and  feasibility  of  this  method 
merit  further  study. 
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3.  Solving  For  The  Eigenvalues 
3.1  The  Eigenvalue  Equation 


Before  the  eigenvalue  equation  can  be  introduced,  the  eigenfunction 
differential  equation  must  be  rewritten  in  more  suitable  form.  To  do  so,  let 
Y$(x)  satisfy 


jo2(x)r,”(»)  +  M*)i7(*)  +  =  o 


Multiplying  (3.1)  by  p(x)  gives 


where  p(x)  and  tt(i)  are  given  by  (1.4)  and  (1.5)  respectively. 


Now  suppose  that  Yg  satisfies  the  boundary  conditions 
17(0)  =  0  and  7,(0)  =  1 
Then  define  ol{9)  by 


-  W 


The  eigenvalues  of  equation  (3.1)  are  none  other  than  the  zeroes  of  oj  ,  see  [2, 
Chap  8]  for  further  details. 


3.2  Standardized  Eigenvalue  Problem 


It  is  also  possible  to  transform  (3.1)  into  more  standard  form  using  the 
following  transformation  scheme. 


setting  *  =  Z(x)  =  /  -  u  reduces  (3.1)  to 
o  (<r(«)/2)l/- 


—  +  3{z)^~  +6Y  =  0 

v  “  /V  <y 


where  0{z)  =  {p{z)-~cr\z)} /{cr{z)/ 2)I/2 

4 


r 

f  —~&{v)du 

Putting  Y(z)  —  g(z)y(z)  ,  where  g(z )  =  e° 


gives 


2  +(#-  <i(z))y  =  0 


mmmwmsm 


where  q(z)  —  —$*{2)  +  and  with  boundary  conditions  y'(0)  =  0  and 

y(b)  —  0  with  b  =  Z(r) . 

The  eigenvalues  of  (3.5)  are  the  same  eigenvalues  as  those  of  (3.4)  and  (3.2)  . 
Moreover  the  eigenfunctions  of  (3.5)  are  easily  transformed  into  those  of  (3.2)  . 
In  particular  if  il>n(x)  is  the  eigenvalue  corresponding  to  an  for  (3.5)  ,  then 
<j>n(x)  =  g(Z(x)mZ(x)). 

The  following  powerful  asymptotic  results  (as  n  — *■  00  )  are  known  about  the 
eigenvalues  and  eigenfunctions  of  the  standardized  problem  (cf  (10,  p.  19]). 


an  =  nV/'i2  +  0(1) 

(3.6) 

ib„(x)  =  (2/b)1/'2cos(n7Vx/b)  +  0(— ) 

n 

(3.7) 

i/Sn(x)  =  —n(2/b3)1^2sin(mcx/b)  +0(1) 

(3.8) 
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4.  Obtaining  First  Passage  Moments 

To  solve  (2.1)  we  need  to  produce  the  moment-sequence 
{A Ik{x,r)  ,  1  <  k  <  rc},  three  such  methods  are  outlined  below. 

4.1  Complex  Integration  To  Invert  Moment  Generating  Function 

Corollary  6.4  shows  that  ^g(x,r)  is  an  analytic  function  of  9  around  0,  and 
so  Cauchy’s  formula  implies  the  identity 


K{x,r)  i  v*x,r) 

n!  27T»'  |tf| lg0  9n+ 1 

where  90  is  sufficiently  small. 

The  integrals  in  (4.1)  may  be  computed  numerically  using  Gaussian  quadrature 
to  minimize  the  number  of  values  of  9  to  be  evaluated,  and  then  taking  the  real 
part.  This  reduces  evaluating  the  equation  to  calculating  a  small  number  of 
values  of  ^(z,r)  .  Evaluating  %(x,r)  may  be  done  by  either  using  finite 
differences  to  solve  the  boundary  value  problem  (cf  equation  (7.6)),  or  by  using 
the  series  method  suggested  in  the  differentiation  approach,  or  some  hybrid  of 
series  and  finite  differences.  An  important  virtue  of  estimates  of  .V/„(z,r)  based 
on  formula  (4.1)  is  that  the  accuracy  of  these  estimates  is  independent  of  the 
accuracy  of  the  n — 1  smaller  moments,  unlike  the  methods  given  in  sections 
(4.2)  and  (4.3)  below. 

4.2  Recursive  Integration 

This  approach  iteratively  uses  (6.1)  to  compute  successive  moments.  This  is 
feasible  when  the  successive  moments  form  a  closed  family  of  integrals  (compare 
example  1),  or  when  only  a  few  moments  are  desired. 

Example  1 

Choosing  parameters  <r(x)  =  2{x  +  <70)  and  /i(j)  =  v  where  u  #  0 
gives  rise  to  the  iteration: 

T  Ul 

Mn{x)  =  f  -  —  —  J  Mn_l(u)(u  +  croy~ldudvj 

z  [  W  OqJ  o 

An  easy  induction  will  show  that  for  u  not  integer  \ln_x(x)  satisfies  the 
expansion 

K-i(z)  -  c[»,0]  +  £  c[n,£](z+a)*  4-  £  rf[n,*j(z+a  )*“" 


%(x,r) 


WW 


where  the  coefficients  are  calculated  using  the  iteration: 

c  [n,£J  =  —nc[n—l,k—l\/(k(k—l)+vk)  where  k  >  1 

d\n,k]  =  — n  </  [  n  — 1 ,  £  — T  ]  /(( —v+k )( — —1  )+*>( —v+k ) )  where  k  >  2 

Finally  c[n, Oj  and  d[n,l]  are  determined  by  solving  the  two-dimensional  linear 
system  arising  from  the  boundary  conditions  Mn(r )  =  0  and  Mn'( 0)  =  0  . 

For  integer  v  ,  closed  formulas  for  all  moments  are  obtainable,  but  the 
calculations  will  be  messier. 

4.3  Differentiating  the  Moment  Generating  Function 

Successive  moments  may  be  obtained  by  calculating  the  0-derivatives  of  the 
moment  generating  function  'I '$(x,r)  at  6  —  0.  This  approach  is  facilitated  by 
Kent’s  observation  in  [7]  that  4>$(x,r)  =  Y9(x)/Tg(r)  where  T^x)  satisfies 

^(z)T,"(x)  +  /i(x)T/(x)  +  6Te(x)  =  0  (**-2) 

with  initial  conditions 

T,'(0)«0  ,  T,(0)#0  (4.3) 

To  solve  for  the  0-derivatives  of  ^(x.r)  it  suffices  to  solve  for  the  derivatives 
of  T*(x).  We  can  obtain  Ttf(x)  using  the  series  expansion  method  around  x=0  to 
solve  (4.2)  .  Under  certain  regularity  conditions,  the  Taylor  series  coefficients 
may  be  differentiated  with  respect  to  6  ,  and  the  series  summed.  This  process  is 
illustrated  in  Example  2  below. 

Example  2 

We  indicate  how  the  technique  in  (4.3)  may  be  applied  to  Example  1: 


T,(x)  =  £  hj(9)xJ 


j-  o 


Equation  (4.2)  implies  that 


cr(x) 

7-0 


7+2 

o 


b^,{6)x]  +  M-r)  £  (7  +l)0J  +  i(0)x;  4-  0  £  bj(9)x>  =  0 
7-o  7-o 


Since  n{x)  =  u  and  o(x)  =  2(x+<r0)  we  deduce  that 

9bj{9)  +  (v+j(j+l))bj+l(9)  +  (T0(j  +2)(j  +l)b}+2(6)  =0  for  j  >  1 


(4.4) 


where  the  boundary  conditions  are 
bo(0 )  =  1  and  b^d)  —  0 


Repeatedly  differentiating  (4.4)  will  give  successive  iterative  formulas  for 
computing  T^n\x)  .  For  example  for  n=l 

9b/{9)  +  bj(6)  +  (u  +  j(j+l))bJ+l'(d)  +  cr0(j+2)(j+l)bj+2'(9) 

With  initial  conditions  b0'(9)  =  6/(0)  —  0  . 

If  the  above  iteration  diverges,  we  can  always  try  renormalizing  by  x  and 
calculating  bj'\9)xn .  If  r  is  sufficiently  small  then  renormalization  will  suffice, 
otherwise  T^(z)  can  be  calculated  by  successively  moving  out  from  0  towards  r 
as  suggested  in  section  (5.1)  below,  and  then  using  the  contour  integration 
method  suggested  in  section  (4.1). 


5.  Some  Remarks  About  Computing  ^(x.r) 


5.1  Computing  ^(O.r) 


The  series  expansion  method  may  not  permit  solving  for  4^(0 ,r)  in  a  single 
step.  However,  suppose  the  series  converges  for  some  y  e  (0,r)  ,  ie,  it  is  possible 
to  compute  by  the  series  method.  We  may  use  4^(0,  y)  as  a  bootstrap 

to  calculate  4^(y ,r)  as  follows.  Observe  that  4^(0, r)  =  4^(0,y  )4^(y ,r)  .  Thus 

fy(o,y)  d4^(y,r) 

5^(0.  y)  dy 


dy 

Using  this  initial  condition  and  Kent’s  normalization  technique,  it  is  possible  to 
calculate  4^(y ,r)  starting  from  y  rather  than  from  0. 

5.2  Interpolating  4^(x,r)  And  A  Related  Boundary  Value  Problem 

Suppose  that  ^(x,r)  and  4^(y,r)  have  been  obtained  (x  <  y)  ,  and  it  is 
desired  to  calculate  4^(,z,r)  for  z  6  ( x,y )  .  The  multiplicative  character  of 
4^(x,r)  implies  that 


%(z,r)  -  %(z,y)%(y,r) 


It  thus  suffices  to  determine  V9(z,y)  .  We  have  %(x.y )  =  4'tf(x,r)/4'*(y,r)  and 
=  1  •  Therefore  it  suffices  to  find  h9(z )  (  =  ^(^,y))  such  that 

j<r{z)h9"{z)  +  n(z)h9’{z)  +  6h„(z)  «  0  (5.1) 


with  boundary  conditions  /i*(x)  =  4'5(x,r)/4'„(y,r)  and  h9(y)  =  1  .  These 
boundary  conditions  uniquely  determine  h9  .  To  solve  for  h9  .  first  find  £0  and 
satisfying  (5.1)  .  where  £,'(x)  =  1  -i  and  £,(x)  =  i  ,  for  i  =  0.  1  .  Then  set 


h9(:)  =  %(x,y)^(z)  +  Uz)<l 


kCl 


(*))/&(*) 


8.  Theoretical  Complements 


In  this  section  some  of  the  properties  of  moments  of  rz  are  examined,  but 
first  some  new  notation  is  introduced. 


Define  Mn(x,y )  =  Ez(t *)  for  x  <  y  and  n  >  0  ,  and  let 
dMn(x,y) 


Mn'(X’y)  1 

Mn"(x,y) 


dx 

d2Mn{x,y) 


dx 2 


6.1  Recursive  Equations  For  Moments  Of  Tx 


The  functions  Mn(x,y)  jointly  satisfy  the  iterative  differential  equation  (cf 
[5],  p.  203,  equation  (3.38)  ). 


1  Jlr 


j<r(x)Mn"(x,y)  +  M(x)M„'(x,y)  +  nMn_t(x,y)  -  0 


(6.1) 


subject  to  Mn'(0,y )  =  0  and  Mn(y,y )  *  0  . 


6.2  Lipschitz  Conditions  For  Moments  Of  Tz 


Lemma 

Mn(x,y)  is  a  smooth  function  in  x  and  y  jointly  ,  and  there  exists  a  constant 
C  such  that 


|A/n(x,y)l  <  Cn  y2n~\y-x)n\ 
and 

|.V//(j,y)|  <  Cnr‘n-1)n! 
for  all  x  <  y  . 

Proof: 

Set 


(6.2) 


!{x)  -  exp{-f^&  d£  } 
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and 


m(x)  *=  l/[<r(x)s(x)]  . 
Rewriting  (6.1)  as  follows 


d  Mn'(x,y ) 
dx  s(x) 

implies  that 


=  -2nMn_l(x,y)m{x) 


y  v 

Mn{x,y )  -  2 nj  [/Mn_1(^,y)m(^)d^5(77)d7 7 

x  0 

By  virtue  of  continuity  there  exists  a  constant  K  such  that 
ll«  II  =  Sup  k(x)|  <  K 

0  <  x  <  r 

and 

Ik  II  <  K  • 

Therefore 


y  v 


|V„(i,»)|  <  2nJ  f  K*II.V„_,lld(dn 

z  0 


=  2 n  /t:2||iV/n_l||y(y-x)  <  2 K2y2n  ||A/„_i|| 


An  easy  induction  implies  that  Mn(x,y )  is  a  smooth  function  in  x  and  y  . 
moreover 

l|.v„(*,*)||  <  2' K-"  !!■"-'{, J-X)n\  (6.3) 

and 

||M/(x)2/)||<2VC2V(”-1)n! 

Taking  C  =*  (2AT2)  completes  the  proof. 

(6.4)  Corollary 

%(x,y)  <  00  whenever  \6\  <  C~ly  2,  and 

-  E  ~ 7  Mn(x,y ) 
n-0  n! 


(6.4) 


6.3  Infinitesimal  Relations  Governing  First  Passage  Moments 


Proposition 

Define 


Un{x,y) 


dMn(x,y) 

dy 


and 

«»(*)  33  Un (x,x) 

Then  M„(x,y )  ,  Un{x,y)  ,  un(x)  satisfy 


K{o,y) 


n 

j 


Mj(0,x)Mn_j(x,y) 


(6.5) 


j-o 


K(  0,x)un_j{x) 


(6.6) 


Proof 

Conditional  on  Jf(0)  =  0  the  strong  Markov  property  (SMP)  implies  that 
{  Tal  >  Ta,  —  ral  >  ~  ra.  }<  where  0  <  at  <  a2  <  •  •  •  <  a,  <  r,  form  a  set  of 

independent  random  variables.  In  particular 


m  =  £®[(t,  +  r,  -r,)*) 


j 


E° |  r,)>  1  £*[(r,-r,)»->  ] 


which  coincides  with  (6.5)  .  Using  a  litcle  algebraic  manipulation  on  (6.5)  shows 
that 


n  —  1 


[Mn(0,y)-M„(0,a;)]/(y-i)  -  £ 


j-o 


Mj{0,x)Mn_j{x,y)/{y-x) 


(6.7) 


Now  letting  y—*x  in  (6.7)  yields  (6.6)  .  QED 
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Comments 

Equations  (6.5)  and  (6.6)  provide  son:  nice  intuition  about  the  way  that  first 
passage  times  from  x  to  y  depend  on  fi.  passage  times  from  0  to  x. 

Equations  (6.5)  and  (6.6)  are  similar  to  (6.1)  ,  but  may  capture  better  the 
dependence  of  higher  moments  on  lower  moments.  From  a  practical  point  of 
view  equation  (6.1)  is  certainly  preferable  for  moment  calculation.  In  section  4, 
other  methods  are  proposed  for  calculating  the  moments  Mn(x,y)  . 


7.  A  Representation  Result 


Theorem 


V$(x,y)  «  exp{  £  /  «»(*) 

n— 1  r  '*• 

where  {un(x)  ,  n>  1}  satisfy 


yoVWOO  +  mOO^iOO  -  1  . 


and  for  n  >2 


^<r(x)un\x)  +  -  —o2^)  E 

4  *-i 


*!(«-*)! 


Proof 

We  begin  by  showing  that  (7.1)  holds  for  ^#(0,r)  . 

Let  xj  »  jr  /L  for  0  <  j  <  L.  Using  the  SMP  as  in  the  proof  of  (6.5), 

*#(0,r)  -  n1  %(xi^ l) 

j-  o 

Taking  logarithms, 

logfy(0,r)  -  £  l°g  %(xj,xJ+l) 
j-o 


Applying  proposition  (8.4)  to  the  above  yields 

log^(O.r)  =  £  Ify(*/»*/+i)  “  1  +  ( )-lfg{ ^(x; ,x>+ 1 ))]  . 

j-o 

Since  ^(x,!)  =  1  it  follows 

ty$(xj ,jy+1)  -  1  =  -^—(xj,Xj  +  aL  x)  L  1 

where  0  <  a  <  1  .  Also  since  tyt(x,y)  is  a  smooth  function  in  x  and  y 
it  follows  that  — —  ’  ■  is  uniformly  bounded  on  the  region  0  <  x  <  y 
there  exists  a  function  A(L)  such  that  A(Z»)~0(1) 

log^(0,r)  -  *£  (*#(*/, *y+1)-l]  +  ^(L)L'1  . 
j-o 


(7.1) 

(7.2) 

(7.3) 


jointly  , 
<  r.  So 


Replacing  each  9f(j.y,xy+1)  in  the  above  by  the  expression  in  (6.4),  yields 


log'I'f(C) 


L- 1  oo  0niVfn(z,,x,+1) 

.<■)-  s  E  — *  ;  +  0(t-‘) 

7-0  n-1  n! 


Since  the  summands  are  positive  the  order  of  summation  may  be  permuted  to 
get 


,r)-  E  E  - -  „  '  +  0(L-') 


log^(0 

n  — 1  /— 0 

The  inner  summation  can  be  expressed  as  a  Riemann  sum 
oo  L-i  A/_  i 

iog^(o,r)  -  s  0n  E  ~  ,  r  -i  T  +  °(L  )  (7.4) 

n— i  y-0  n!  L  ^ 

Equation  (6.2)  implies  that 

^n(xVX)+l)  <-  y-,n/  l2"  <  C  n 

n!  L“1  “  ■ 

where  C0  ■■  Cr2  . 

The  Lebesgue  dominated  convergence  theorem  (applied  to  the  product  space 
{l,2,..}  X  [0,r]  endowed  with  product  of  the  counting  measure  with  the 
Lebesgue  measure  on  [0,r])  implies  the  right  side  of  (7.4)  approaches  the  limit 

log’MO.r)  =  S  ~f  «.(*)* 

n-1  0 

as  L  — ►  oo  ,  or 

^(0,r)  =  exp{£  ~f  un{z)dz} 

n-1  n-  0 

Due  to  the  multiplicative  nature  of  '^>$(xfy)  it  is  easy  to  show  that 

%(x,y)  -exp{£  /  un{3)  — dz }  (7.5) 

n-l  x  n • 


The  function  V$(x,y)  satisfies  the  following  differential  equation  (cf  [5]  pp  203  ) 

(7.6) 


1  _2/  s  52^,(x,y)  #  %  dVt{x,y) 

2^{x)  a?  +*{x)~x  - 0 
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*'  y  -t  ,y. 


dVt(0,  y) 

subject  to  - - - «■  0  and  ^(y,y)  —  1  . 

ox 


Substituting  equation  (7.5)  in  (7.6)  ,  shows  that  the  exponent  in  (7.5)  satisfies 
the  differential  equation 

t 7“»0O)2  -  E  T7«,'(*)]  -  M*)E  t7»„(*)]  +  0  -  °- 


n-1  n! 


«-l 


n  —  1 


Since  the  above  series  converge  absolutely  for  9  sufficiently  small,  we  can 
rearrange  terms  to  obtain  a  single  power  series  in  9  .  Since  this  power  series  is 
zero  for  9  sufficiently  small,  all  its  coefficients  must  be  zero.  Setting  the 
coefficients  of  9  to  zero  yields  equations  (7.2)  and  (7.3)  . 

The  initial  condition  that  °’v -  *  o  in  (7.6)  implies 


tin(0)  =  0  ,  where  n  >  1  . 


R 


8.  Appendix 

8.1  Spectral  Representations  For  First  Passage  Time  Distributions 

Theorem  The  right-hand-side  of  (1.3)  is  the  unique  function  satisfying  (1.2), 

duu  i  oc  1 1 

jointly  continuous  in  x  and  t  on  [0,r]  X  (0,oo)  with  - absolutely 

at 

integrable  over  [0,  r ]  X  (N~l,  N)  for  all  N  >  1. 

Proof: 

Suppose  that  w(x,t)  satisfies  equation  (1.2)  and  the  integrability  conditions. 
Observe  that  for  fixed  t  >  0  w(x,t)  is  a  continuous  function  of  x  belonging  to 
L'(p)  .  Therefore  w(x,t )  possesses  the  orthogonal  expansion 

w(x,t)  =  £  c*(0  Mx) 

k-l 

where 

r 

ck{t)  -  /  w(x,t)  <t>k{x)p{x)dx 
o 

Multiply  (1.2)  by  <f>k{x)  and  integrate  over  [0,r]  to  get 

/  </>k(x)p(x)dx  =  /  Aw(x,t)<i>k(x )p(x)dx  . 

o  ol  o 

Now  since  A/(x)p(x)  =  tt(x)/"(x)  +  ^{x)J\x)  (cf  (3.1)  and  (3.2)),  a  simple 
integration  by  parts  shows 

r  r 

f  Aw(x,t)d)k(x)p(x)dx  =  J  w(x,t)A  4>k{x)p{x)dx 
o  o 

The  relationship  A  dk(x)  =  —akdk(x)  implies  that 

^  3  w  ( x  t }  r 

I  <t>k{x)p{x)dx  =  -  Jakw(x,t)dk(x)p(x)dx 

o  ac  o 

Integrating  both  sides  with  respect  to  t  over  [u0,«]  and  permuting  the  order  of 
integration  on  the  left-hand-side  yields  (permissible  because  Fubini's  Theorem 
applies  to  the  absolutely  integrand  ^  cf>k(x)p(x)). 

r  u  dwix  t)  u  r 

I  f  — “-^(x^x)  dt  dx  =  /  /  -akw{x ,t)4>k{x)p{x)  dx  dt 

0  «o  u0  0 


Using  the  definition  of  c*(£)  on  the  right-hand-side  of  the  last  equation  implies 

u 

ck(u)  +  C  —  /  ck(t)dt  , 
u0 

where  C  is  an  arbitrary  constant. 

Therefore  ck(t)  «=  cke~°‘it .  It  remains  to  determine  the  constants  ck  ,  but  they 
may  be  derived  from  the  boundary  condition  w(x,0)  =  1  for  0  <  x  <  r  as 
follows.  Since  p  is  continuous,  w(x,0)  =  1  for  almost  all  p(dx),  and  1  €  L~(p). 
So 

1  -  S  ckM*)  M 

*-i 

where 


r 

ck  —  J  <t>k(x)p(x)dx  (8.2) 

o 

Now  Theorem  1.9  of  [10]  implies  that  the  right-hand-side  of  (8.1)  converges 
pointwise  to  1  on  (0,r)  .  Therefore  w{x,t)  has  the  representation 

w{x,t)  -  £  ck  e~aktd>k(x)  (8.3) 

k- 1 


To  prove  the  converse,  suppose  that  w(x,t )  is  defined  by  (8.2)  and  (8.3)  jointly. 
Equations  (3.7)  and  (8.2)  imply  that  the  coefficients  ck  ,  k>  1  are  uniformly 
bounded.  Hence  for  t  >  e  >  0  the  series  converges  uniformly  to  a  function 
continuous  on  the  product  [0,r]  X  (e,oo)  .  The  uniform  convergence  and 
boundary  conditions  on  the  eigenfunctions  imply  the  boundary  conditions  on  x 

.  The  integrability  conditions  on 


follow  in  similar  fashion.  The 


boundary  condition  w(x.O)  =  1  for  x  E  (O.r)  follows  from  Theorem  1.9  of  [10] 
and  equation  (8.1)  The  differential  equation  (1.2)  may  be  derived  from  the 
definition  of  derivatives  as  limits  of  divided  differences,  and  the  dominated 
convergence  theorem  applied  to  series.  QED. 


It  should  be  noted  that  the  first  passage  time  distribution  satisfies  the  regularity 
conditions  of  the  theorem,  and  therefore  must  have  representation  (1.2). 


8.2  Convergence  Of  The  Finite  Approximations  To  The  Infinite  Vector 

It  will  now  be  shown  that  the  solution  vector  to  system  (2.1)  ,  denoted  by 


P**1  **  {p*(,) >  k  <  k  <  n},  converges  at  rate  n  2  component-wise  to  the  infinite 
vector  p  =  {pk  ,  k  >1}  ,  where  p*  =  c*<^(;r)  . 

Theorem 

Ip”  -  p. |  -  0(A) 

n~ 

Proof:  Define  v*"1  ,  Slm)  and  rji']  as  follows: 

n 

i^k')  =  £  VjPj  ^here  0  <  k  <  n— 1 
00 

—  mt  —  ^*’  =  5]  fijpj  where  0  <k  <  n-1  (8.4) 

J-n  +  l 


=  pt  -  pt«*» 

Observe  that  {u^  ,  1  <  k  <  n}  is  the  solution  to  (2.2)  where  {mn  ,  1  <  it  <  n} 
has  been  replaced  by  { 6 i*1  ,  1  <  k  <  n}.  In  particular 

l>i  =  £  hj  gn—j{tii,P2,  •  ■  ,P*-i  ,  Mt+i,  •  •  p„)  /  (My  —  A*/t )  (8.5) 

l<J<n  1  <  j  <  n 

;>r 


For  1  <  j  <  n 

'  •  '  ./'t-i  ,  P*+i,  •  •  ■  ^j)l  <  Hk~J~k 

•  Also  equations  (3.2),  (3.6)  and  (3.7)  together  imply  that  |cn  |  =  0{— )  ,  thus  (8.4) 
and  (3.6)  imply 

|<5;i  <  Cn~‘-> 


>o 


Is  Vn-jiiv  1^2, 
1 


;-i 


<  »rk  £  -  vrk-l<K- T) 

/— l  n* 


Thus 

I  £  9n— ;((P  1  >A*2» 


.P*_,,P*+l)  ••  -p,))5/”|- 0(^-7— ) 


ft—  /t  —  1 
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Finally,  the  denominator  of  (8.5)  may  be  written  as 


n  -»k)  =  Hk~k4n) 

1  <  j  <  n 

j+ir 

where 


4n)-  n  (<*/-/*»)  n  a -•£*-) 

1  <  j  <  k- 1  /fc  +  1  <  j  <  n  rk 

Now  equation  (3.6)  and  the  Weierstrass  Product  Convergence  test  jointly  imply 
that 

lim  =  d  >  0 

n—+oo 

Now  it  follows  that  rj**’  =*  0(-“)  as  claimed. 

n~ 


8.3  A  Conjecture 

It  is  interesting  to  note  that  pk[x)  is  linearly  proportional  to  the 
eigenfunction  <f>k(x).  and  therefore  Apk(x)  =  akpk(x).  This  suggests  that  pk''(x) 
will  approximately  satisfy  this  relation.  Observe  that  Amk(x)  —  ~  mk_l(x). 
Thus 

Apk'\x)  =  £  •  •  •  ,/<*-,  ,  pk+l,  ■  ■  •  pk)  /  II  (Vj  -  Vk 

1  <  j  <  n  —  1  •  1  <  j  <  a 


Comparing  this  with  (2.3)  suggests  that 

gn-i-i(Uli<fu.- '  •  -  •■*<„))  i 

lim  - — : -  =  c*k  =  — 

1— a„_;((M |.M2,  -Mk-i  M,))  h;- 


8.4  Proposition 
For  1  <  i  <2 

|/o«7(z)  —  x  -f  1 1  =  (x— l)‘g(x)/2  ,  where  |p(x)|  <  1 

Proof:  The  mean  value  theorem  applied  to  log(x)  at  x=l  gives: 
log{x)  *  log(\+x—\)  =log(  1)  +  (x— 1)  —  (x— l)~/2(l+a(x— 1)~2) 

where  0  <  a  <  1  .  The  prop  now  follows  from  x  >  1  and  a  >  0  . 
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